Influence of non-universal effects on dynamical scaling in driven polymer translocation 
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We study the dynamics of driven polymer translocation using both molecular dynamics (MD) 
simulations and a theoretical model based on the non-equilibrium tension propagation on the cis side 
subchain. We present theoretical and numerical evidence that the non-universal behavior observed 
in experiments and simulations are due to finite chain length effects that persist well beyond the 
relevant experimental and simulation regimes. In particular, we consider the influence of the pore- 
polymer interactions and show that they give a major contribution to the non-universal effects. In 
addition, we present comparisons between the theory and MD simulations for several quantities, 
showing extremely good agreement in the relevant parameter regimes. Finally, we discuss the 
potential limitations of the present theories. 



I. INTRODUCTION 

The translocation of a polymer through a nano-sized pore occurs as a part of many biological processes and functions, 
such as DNA and RNA translocation through nuclear pores, protein transport across membrane channels and virus 
injection pQ. The translocation process is also envisioned to have several biotechnological applications, including rapid 
DNA sequencing, gene therapy and controlled drug delivery [2J. Recently, the hope of realizing a third generation DNA 
sequencing method using nanopore sequencing devices has prompted rapid advancement in experimental methods and 
in technological applications 0] , leading to commercialization of the nanopore sequencing technology in the near 
future [5]. 

However, despite the technological advances and considerable experimental [SJIHIIT] and theoretical [5 HT71 H^l 1201 1221 - 
1251 |3"U1 15T1 154TI31)] efforts to understand the basic process, until very recently the fundamental physics of driven polymer 
translocation has remained elusive. The first attempts to model polymer translocation theoretically were made by 
Sung and Park [9] and later by Muthukumar [10], who viewed polymer translocation as a one-dimensional barrier 
crossing problem of the translocation coordinate s (the length of the subchain on the trans side) , with the activation 
barrier determined by the free energy of the polymer derived from equilibrium enumeration of random walks. However, 
if the process is facilitated by an external bias, the process is inherently out of equilibrium due to the long relaxation 
time of the polymer chain. This fact was first recognized by Kantor and Kardar |12j , and was later studied extensively 
in simulations [T2J [271 HFJ EDI 131], with chain lengths and driving forces in the typical experimental regime. It was 
observed that with increasing driving force, significant non-equilibrium effects come into play |301 134] . The driven 
translocation problem was also studied by several authors with different approaches, such as the fractional Fokker- 
Planck method by Metzler and Klafter [11] , scaling theory supplemented with the fractional Fokker-Planck method by 
Dubbeldam et al. [13] and by Vocks et al. using a method based on local memory effects near the pore [2]. However, 
the first truly non-equilibrium treatment was given by Sakaue, using an ingenious tension propagation theory based 
on force balance, mass conservation and self-similarity of the polymer [15HT7| IT9]. However, the theory was solved in 
the asymptotic limit of long chains, neglecting, e.g., the interaction between the pore and the polymer. It was pointed 
out already by Storm et al. [7] that this interaction may significantly alter the dynamics of driven translocation, which 
was further demonstrated in simulations [29, 32, 33, 45 . Hence, the agreement of the theories with experiments and 
simulations remained modest. 

Recently [43] , we have generalized Sakaue's tension propagation (TP) model of driven translocation for finite chain 
lengths and included the pore-polymer interaction as an additional friction force. It was shown that the model is 
in excellent agreement with high-accuracy molecular dynamics simulations. In addition, we showed that due to the 
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pore friction and finite size corrections to the tension propagation equations, the asymptotic limit is well beyond the 
present computational capabilities and simulation algorithms. For example, the scaling exponent a, which relates 
the mean translocation time t to the chain length N as r ~ Nq, actually retains a fairly strong dependence on the 
chain length even up to Nq sa 10 5 . Therefore, the scatter (non-universality) of a, widely studied and reported in the 
literature, is in fact a finite chain length effect. 

The purpose of this work is to study further the non-universalities in driven polymer translocation brought on by 
finite chain length effects, and to see how the asymptotic limit is approached in different regimes. To this end, we 
use a combination of theoretical methods and molecular dynamics (MD) simulations. For the former, we use the 
theoretical Brownian dynamics tension propagation (BDTP) model introduced in Ref. [33] , which is based on the 
tension propagation description and incorporates the finite size effects, and for the latter, we use full A^o-particle 
Langevin thermostatted MD simulations. The BDTP model and its numerical solution is described in Section [TTJ 
while the details of the MD simulations are included in Appendix A. 

The central assumption of the tension propagation theory is that the effect of the trans side subchain on the 
non-equilibrium dynamics is small, and therefore its contribution to the effective total friction is neglected. This 
approximation is very good in the parameter regime typically used in experiments and simulations. However, in 
certain cases, the influence of the trans side chain on the dynamics may be non-negligible. Therefore, we do not 
expect perfect agreement with the BDTP model and the MD simulations. The limitations of the model are discussed 
in Section |IIf| However, it will be shown that for most relevant situations, the BDTP model gives extremely good 
match with MD. In addition, it allows one to go beyond the chain length regime available to MD simulations (No « I0 3 
for high driving force) , and to see how the finite chain length effects diminish as the asymptotic limit is approached. In 
addition, in Section [Hi] we present detailed comparison between the theory and MD simulations for the time evolution 
of the translocation coordinate s(t), the scaling of translocation time with the driving force and the monomer waiting 
time distribution. In our previous work we have explained and demonstrated the applicability of the BDTP model 
to both two and three dimensional systems |43j . In this work we provide new results for three dimensional geometry, 
some of which are directly relevant to gain better understanding of the experimental results for DNA translocation 
through nanopore. We also relate our findings to previous MD simulations found in the literature, and to the theories 
of Refs. [T5H3Z1Q1M]. 



II. MODEL 



A. General framework 



The purpose of the BDTP model is to present a coarse-grained, minimal model of driven polymer translocation [43j . 
Since the driven translocation is a complex, non-equilibrium dynamical process, the rigorous solution of the full 
problem from first principles seems impossible at present. Instead, the BDTP model presents a phenomenological 
description, which interpolates between the low-force and high-force regimes. In the low-force regime, the fluctuations 
and the chain entropy become important. Hence, as a general framework, we adopt the description used by Sung 
and Park [3] and Muthukumar |10j . the one-dimensional barrier crossing problem of the translocation coordinate 
s. Here, the chain starts from the cis side with one end inside the pore and is considered as translocated once 
s = aNo, with a the segment length. The free-energy due to chain entropy and the chemical potential difference Afj, 
is J"0) = (1 - l')k B T\n [f (N - £)] + |A)ti. Here 7' is the surface exponent (7' = 0.5, « 0.69, « 0.95 for an 
ideal chain, and a self-avoiding chain in 2D and 3D, respectively), and fc^T is the thermal energy. From T(s), the 

Brownian dynamics equation for s in the overdamped limit follows as = (I — 7')fcsT aN > _ s — \ — + ((t)- 

Here r is the effective friction, and ((t) is Gaussian white noise satisfying (£(i)) = and (C(*)C(*0) = 2FfesT5(t — tf). 
In this framework, the non-equilibrium memory effects at larger driving forces are then taken into account by allowing 
the effective friction to depend on time, T = T(t). This time-dependence is then solved from the tension propagation 
(TP) formalism. 

The central idea of the tension propagation theory is to divide the subchain on the cis side into two distinct 
domains [T5HT71 H9], The first domain, closer to the pore, consists of all the monomers that are pulled towards the 
pore by the external driving force. The second domain consists of the remaining monomers, which are at rest (on the 
average). As the driving force is applied at the pore, the chain begins to move in stages, with the segments closest to 
the pore being set into motion first. A close analogue is a coil of rope pulled from one end, which first has to uncoil 
and become tense before starting to move as a whole. To keep track of the moving part of the chain, one defines a 
tension front, which divides the chain into the moving and nonmoving domains. The front propagates in time as parts 
of the chain further away from the pore are set in motion, as shown in Fig. [l] As the monomers enter the front and 
start to move, the effective friction T increases due to the increased drag between the polymer and the solvent. After 
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FIG. 1: (Color online) The time evolution of the polymer configuration during the translocation process, with time advancing 
from left to right. The arc denotes the position of the tension front R, which separates the chain into the moving and nonmoving 
domains. The last monomer inside the front is denoted by N and the number of translocated monomers by s. 



a certain tension propagation time, t tp , the front reaches the end of the chain and the tension propagation process 
stops. After this time, the chain as a whole is pulled towards the pore. During this stage, the overall length of the 
subchain on the cis side decreases, which reduces the effective friction T. This stage continues until the last monomer 
reaches the pore and finally translocates at time t = f. 



B. Coarse-grained equations of motion 



In the BDTP model, the translocation process is described using only two degrees of freedom. The first one, the 
translocation coordinate s, describes the chain's position with respect to the pore, while the second one, the location 
of the tension front R, describes the response of the cis side chain to the driving force. To both of these degrees of 
freedom one has a corresponding equation of motion, which are coupled via the effective friction T. In the extremely 
high driving force limit, the location of the tension front can be described by only one equation, while in the more 
general case considered in this work, one needs to solve a system of several equations. The derivation of the BDTP 
model is presented in Ref. [43] and the tension propagation theory in Refs. [15HT7]. IT9] . Therefore, in this work we do 
not reiterate the full derivation of the model, but rather present the resulting equations and outline the method of 
their numerical solution. 

For convenience, we first introduce dimensionless units denoted by the tilde symbol as X = X/X u , with the unit of 
length a u = a, force f u = hsT/a, time t u = rja 2 /kBT, velocity v u = a/t u and friction rj u = r), where r] is the solvent 
friction per monomer. In these units, the Brownian dynamics equation for s = s/a reads 



f(i)§ = (i-V) 



1 1 



N -s s 



+ ;+c(t) = ; tot , (i) 



where T = T/rj and, for simplicity, we have assumed that the pore length l p — a. Generalization of Eq. to 
different pore lengths is straightforward (see, e.g., Ref. [21] for a similar case). The dynamics of Eq. is essentially 
determined by the effective friction T, which therefore must comprise all the dominant contributions of the non- 
equilibrium dynamics of the full polymer chain. Although under sufficiently large /, also the trans side subchain 
assumes a highly out-of-equilibrium configuration, it has been shown in Ref. [43] . that for most cases the influence 



of the trans side chain on the dynamics is extremely small. This matter will also be further elucidated in Section III 
of this work. Therefore, to a very good approximation it is sufficient to consider the drag force due to the cis 
side subchain and the frictional interaction between the pore and the polymer. Although the latter is negligible for 
infinitely long chains, it can significantly change the dynamics for finite chains. Formally, we can write T as the sum of 
the cis side subchain and pore frictions, T(t) = fj c i s (t) +fj p . The first contribution can be calculated from the tension 
propagation formalism, while the pore friction fj p has to be determined empirically. In the special case of extremely 
large driving force, one can find ?7 c i s explicitly from the TP equations [43] . More generally, however, it is easier to 
derive the velocity of the monomers at the pore entrance, vq. In such a case, the effective friction is defined as 

f (t) = . (2) 

where do is the line density of monomers near the pore and aoio = ds/dt is the flux of monomers through the pore 
entrance. In either case, determining T(t) essentially reduces to calculating the number of moving monomers, whose 
combined drag force then constitutes the time-dependent part of the friction. 

For moderate driving forces (to be defined rigorously below) the chain on the cis side assumes a trumpet-like shape 
due to the higher stretching close to the point of action of the driving force. In this so-called trumpet (TR) regime, 
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the line density near the pore is <xo = fo , where fo ~ /tot — fjpVo is the force at the pore entrance |39j . To solve 
the velocity vq, one needs to solve a system of four equations: 



dR(t) . 
-g- = v R (t) 



l -A- 1 ' l 'd R {t)- 1 R{t) 1 l 1 '- 1 



cosh b 



R(t) 



(3) 
(4) 



tanh (bd R {tyl^-^ /R) 
**® = tanh(h) ^ (S) 

v (t)R(t)C b = vfl'\ (6) 

where is the Flory exponent, A v relates the chain end-to-end distance to its contour length as R ee — A V ~N V . From 
MD simulations we have measured A v ss 1.15 ± 0.03 in three dimensions. The coefficient Ch = h l^ os H^ with b a 
dimcnsionless parameter that describes the velocity profile on the cis side and ensures global conservation of mass [13] . 
The line density and velocity at the tension front are given by a R and v R , respectively. The numerical solution of 
Eqs. ([3])-@ is described in Section |nC| 



Equations (|3]H(|6| are used when the line density of the monomers at the pore entrance is greater than unity, i.e., 
for fo < 1. For slightly higher forces, the chain assumes a shape consisting of a fully elongated stem followed by a 
trumpet-shaped flower. In this stem-flower (SF) regime, the velocity at the pore entrance is given by 

v (t)R(t)C b = f + v - 1. (7) 

Also in this regime, the line density and the velocity at the boundary are given by Eqs. Q and ([5]), respectively, and 
the time evolution of the front by Eq. ([3]). Note that Eqs. ^ and ^ ensure a smooth cross-over between the TR 
and SF regimes at fo = 1. In practice, we solve Eqs. ([l])-([7]), choosing Eq. ^ over Eq. Q if /o < 1, and vice versa. 

In deriving Eqs. (4|-(7| we have adopted the form v(x,i) = Vq (i) ^^^(b)' 1 ^ ^ or tne ve l° c ity profile of the cis 
side subchain. Here o is a dimensionless parameter that controls the sharpness of the profile and is fixed by enforcing 
global conservation of monomers, i.e., requiring that at the time of translocation s(f) = A^o and R(f) = 0. Although 
the adopted profile is a good approximation to MD simulations (see Appendix B) , the exact functional form of the 
profile is not crucial. The essential requirements are that the profile is non-constant, goes to zero at the tension front 
(x = —R), and allows the global conservation of mass. 

C. Numerical solution of the BDTP equations 

To solve the model numerically we first set the initial values for the translocation coordinate s and the location of 
the tension front R. In the initial configuration, the first chain segment is through the pore entrance and s(i ) = 1. 
This gives the initial condition for Eq. Q. Correspondingly the initial location of the tension front is given by 
R{to) = l//o for fo < 1 and R(io) = 1 otherwise. This gives the initial condition for Eq. 

At the beginning of each time step, the random force is sampled from a Gaussian distribution with the mean 
and variance given in Section |II A| Since the value of s is known at the start of the time step, the value of the total 
force /t t is then also known. The task is then to determine the effective friction T for the next time step. For this, we 
need to find the value for the velocity Vo and the line density do [cf. Eq. p])]. The velocity Vo(t) is given by Eq. ^ in 
the TR regime, and by Eq. (7) in the SF regime. The latter can be solved explicitly for vq, while the former is solved 

numerically by the Newton- Raphson method [40] . Knowing i>o then gives fo — /tot — Vp^o an d = fo • The 
effective friction is then given by Eq. ([2|, and Eq. ([!]) can be advanced in time by the explicit Euler algorithm |42) . 

To solve the effective friction for the subsequent time steps, one also has to find the time evolution of R from 
Eq. ([3i Therefore, one has to know the velocity v R and the line density a R near the tension front. The line density is 
a R = £j[ u 1 , with £h solved from Eq. Q by Ncwton-Raphson iteration. The velocity v R is then given by Eq. |s| and, 
with both v R and a R known, the location of the tension front R for the next time step can be solved from Eq. &3\. 
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We emphasize that the BDTP model is not an alternative formulation of molecular dynamics. It is a model of two 
degrees of freedom, s(t) and R(t). In addition, it has no free parameters that could be used to adjust the results. 
Therefore, the agreement of the model with MD simulations presented in both Ref. [33] and in Section III of this work 
is not self-evident, but a demonstration of the tension propagation formalism's capability to capture the essential 
physics of driven polymer translocation. 



D. Determining the pore friction r\ v 



To complete the BDTP model, we still need to determine the pore friction r\ v . Because r\ p characterizes the 
interactions between the pore and the polymer chain, we do not expect it to have any universal value. Rather, it 
should depend on the geometry and dimensions of the pore and, in a more refined model, on the chemical details 
of the pore and its immediate vicinity. In the context of this work, we only consider coarse grained models of 
the pore, where the pore consists of immobile monomers placed in a configuration that forms the edges of the 
pore (cf. Appendix A). The same approach has been used widely in the literature for MD simulations of polymer 
translocation [2TJ [351 [301 133 [Ml 131] • While it is possible to determine the pore friction also from experimental 
data, in this work we only consider the numerical pores used in our benchmark MD simulations. 

To determine rj p , we look at the movement of the first few monomers and fix rj p in the BDTP model so that the 
velocity matches MD simulations. In this early stage of translocation, where the tension front R is still close to the 
pore, the friction in the system is largely determined by the friction between the pore and the polymer, giving the 
most accurate estimate for r] p . In practice, we look at the waiting time per monomer w(s), defined as the time that 
the individual monomer spends inside the pore. With / sufficiently large, w oc T/f. For small s, the friction T 
and w are mostly determined by f] p . In Ref. [43] . we have measured rj p specifically for the pore geometries used in 
Refs. [23[l34j. for which we had access to the waiting time distributions. In these cases, we have r\ p ~ 5 and rj p w 4 in 
3D and 2D, respectively. Other pore geometries and the effect of r\ p on the translocation dynamics will be discussed 
further in Sec. Ill Finally, it should be noted that for each pore geometry, rj p is fitted only once, as opposed to being 
done separately for each combination of /, r), etc. Thus, r\ p is not a freely adjustable parameter, but is a property of 
the pore. 



III. RESULTS AND DISCUSSION 
A. Waiting time distribution 

We begin the analysis of the results by looking at the waiting time w(s) of individual monomers. We do this by 
solving the BDTP model with parameters typical for molecular dynamics (MD) simulations: / = 5.0, fc^T = 1.2, 
r] = 0.7, v = 0.588 (3D) and pore friction r\ p — 5.0, corresponding to the pore geometry used in, e.g., Ref. [34 . The 
resulting waiting time distribution is an important measure of translocation dynamics. Previously, we have shown 
that the waiting time distribution is reproduced almost exactly by the BDTP model as compared to MD simulations 
in 2D and 3D 43 . In addition, the shape of the waiting time distribution is non-monotonic, with the initial part 
of increasing waiting time w(s) describing the tension propagation stage, and the second part of decreasing w(s) 
being the tail retraction stage (cf. Fig. [2| . It is reassuring to note that the normalized waiting time distribution for 
Ao = 256 looks very similar to that for Ao = 128 for the same parameters reported earlier in Ref. 43 . 

Between the two stages is the maximum of w(s), i.e., the moment of maximum friction T, which occurs when the 
tension front reaches the Ao:th monomer of the chain. For sufficiently large A , the translocation velocity is, according 
to Eqs. ^ and ([7]), Vo(t) cx Immediately after the tension propagation stage, the location of the tension front 

is R oc Nq . Therefore, the maximum waiting time, w max , should scale with chain length as iy max oc [Oo(ftp)] oc Ag . 
This is indeed the case, as shown in Fig. [3j which displays the collapse of the waiting time distributions for different 
Ao onto a single master curve. Since the area under the w(s) curve gives the average translocation time, one has 
r ~ Nq with a « 1 + v. However, even for Ao = 10 5 , the location of u> max slowly moves to the right and the collapse 
to the master curve is not exact. This shows that the chain length A = 10 5 is still not in the asymptotic limit! 
Consequently, for Ao = 10 5 , the scaling of the average translocation time is also not exactly r ~ Ag as we will 
discuss below. 
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FIG. 2: (Color online) Comparison of waiting times w for MD (squares) and the BDTP model (circles). The agreement of the 
BDTP model with MD simulations is excellent, and reveals the two stages of translocation: the tension propagation stage of 
increasing w(s) and the tail retraction stage characterized by decreasing w(s). For both MD and BDTP the data has been 
averaged over 2000 successful translocation events and the system parameters are the same (No = 256, / = 5, ksT = 1.2, 
77 = 0.7). 




FIG. 3: (Color online) The waiting times per monomer w(s) as a function of monomer number s. In the long chain length 
limit, w(s/N ) ~ Nq , implying the scaling r ~ N 1+v . Model parameters used are / = 5.0, k B T = 1.2, rj = 0.7, v = 0.588 (3D). 



B. Dependence of r on the chain length iVo 

Previously the BDTP model has been compared with MD simulations in several regimes, and has been shown to 
reproduce the correct dependence of the average translocation time r on the chain length Nq within the numerical 
accuracy [43] . Some of the results are gathered in Fig.|4j where we show the translocation time exponent a (defined via 
t ~ Nq) for both the BDTP model and MD simulations. The figure shows that merely by using the same numerical 
values for the parameters in the BDTP model and in the MD simulations, the MD results are reproduced by the 
theory with good accuracy. It is also clear that the numerical value of a depends on several parameters such as chain 
length, driving force, friction and pore structure. However, as shown by our analysis of the waiting time distribution 
w(s) and discussed in Ref. [43], all the results shown in Fig. [4] are in the finite chain length regime. Therefore, the 
apparent non-universality of a is not surprising. To further study the finite size effects, we look at the dependence 
of r on the chain length Nq, first for different driving forces / and then for different pore geometries by varying the 
dimensionless pore friction fj p . 
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FIG. 4: (Color online) The exponent a (r ~ No) as a function of driving force / and chain length No. The position of the 
circle denotes the combination of (/, No) and, the size reflects the value of a. Next to the symbols, the numerical values of a 
from the BDTP model are shown in comparison with the values from MD simulations written in parentheses. The asymptotic 
value of a is indicated in the upper right corner. 
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FIG. 5: (Color online) The effective exponent a(No) = j^jj- as a function of chain length Nq for the BDTP model solved 
for driving forces / = 0.75,1.5,2.0,3.0,5.0,10.0. Other parameters are fcsT = 1.2, r\ = 0.7, rj p — 5.0. Errors are of the 
order of the symbol size. Inset: a(No) for / = 3.0,5.0,10.0 up to No = 10 10 , showing the approach to the asymptotic value 
a(No — > oo) = l + f. 

1 . Effect of the driving force 

To study the dependence of the average translocation time r on the chain length Nq, we define the effective (running) 
exponent a(No) = f^Z ■ I n the finite chain length regime, the effective a(No) is a function of the chain length and 
approaches the asymptotic value for long chains. In Fig. pi we show a(No) for different driving forces up to chain 
length N = 10 4 . Throughout this regime, the exponent apVg) shows clear dependence not only on the chain length, 
but also on the driving force /. This is a clear indication of non-equilibrium behavior and finite chain length effects. As 
the chain length is increased further, the different curves approach the same asymptotic value of a(No — > oo) — 1 + v, 
as shown in the inset of Fig. [5] However, the approach is extremely slow: within the numerical accuracy of the BDTP 
model, the asymptotic value is not reached until N « 10 9 . 

Another interesting fact is the dependence of a on the driving force /. Our results show that as / is increased, 
a decreases for any fixed chain length Nq ^ 10 . However, in the literature, there are conflicting reports on the 
dependence of a on /. In Ref. [34) it is reported that a decreases as / is increased, in agreement with our results. 
The BDTP model also agrees with the numerical values of a reported in Ref. [51] with excellent accuracy (see Fig. [4] 
and Ref. [13]). On the other hand, MD simulation studies by Lehtola et al. [5S] and more recently by Dubbeldam 
ct al. [21] report the exactly opposite behavior that a increases with /. However, we argue that neither of these 
studies corresponds to the canonical driven translocation problem, where the chain starts initially on the cis side 




FIG. 6: (Color online) The pore friction rj p as a function of the pore diameter d. The symbols indicate the results obtained 
from MD simulations for different solvent frictions 77, while the solid line is an empirical fitting function. The pore friction r\ v 
includes a non-vanishing contribution pr\ from the p monomers inside and in the immediate vicinity of the pore. For p, we find 
p ~ 2.5, giving the total pore friction r\ v « 2.73/ '(d — 1)?7lj + 2.577, where 77L.1 is given by the choice of the Lennard- Jones units 
(see text). The chain length is iVo = 100, with other parameters the same as in Fig. [3] 
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FIG. 7: (Color online) The effective exponent a(No) averaged over different chain length regimes shown as a function of pore 
diameter d. Model parameters are the same as in Fig. [6] 



and may either successfully translocate or slip back to the cis side, depending on the fluctuations. In Ref. [26] . the 
low-force simulations were performed by placing the polymer chain initially halfway through the pore to facilitate 
successful translocations |46j . Because of the intrinsic non-equilibrium nature of driven translocation, this approach 
does not give results that can be directly compared with theory or experiments, where the process always starts with 
the whole chain initially on the cis side. In Ref. |21j . another method was used to make the low- force simulations 
possible. Here, the authors prevented the chain from escaping back to the cis side by making the first monomer too 
large to fit through the pore. This is equivalent to enforcing an artificial reflecting boundary condition on the first 
monomer that prevents the escape. Although such a boundary condition has been used several times in the literature, 
it fundamentally changes the system's behavior in the low-force limit, as we will discuss below. 



2. Effect of pore size and pore friction rj p 

In the finite chain length regime, it is conceivable that the translocation dynamics is affected by the local neighbor- 
hood of the pore. The effect of the pore size and geometry has been previously studied with MD simulations in, e.g., 
Refs. [29l|4H|45]. It has been shown that the exponent a depends on the details of the pore, although no systematic 
study on the nature of the dependence has been performed. In the BDTP model, the effect of pore geometry is 
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FIG. 8: (Color online) The effective exponent a(No) as a function of chain length No for the BDTP model solved for different 
ratios fj p of pore and solvent friction. Parameters used are the same as in Fig. [6] 

mapped into one parameter, the pore friction rj p . To study the effect of local pore geometry, we have measured r\ p for 
different pore diameters d by running MD simulations (for details, see Appendix A) and matching the waiting time 
distribution with the BDTP model by fixing r\ p separately for each pore diameter d. The results are shown in Fig. [6] 
For the smaller pore, the interactions between the pore and the polymer are stronger, with rj p increasing rapidly as 
the diameter of the pore approaches the size of the monomer, fn addition to the friction between the pore walls 
and the polymer, r\ p includes the effective drag force of the monomers inside the pore and in its immediate vicinity. 
The latter constitutes a non-vanishing contribution to the pore friction, and is extracted from Fig. [6] by collapsing 
the data for different r\ on a single master curve. Empirically, the total pore friction is thus approximately given by 
r/p w 2.73/(g? — + 2.577, where the first term describes the interaction between the pore walls and the polymer, 

and is similar in form to the one used by Storm et al. [TJ. The factor is given by the energy (e) and length scale 
(a) of the Lennard- Jones interaction and the monomer mass m (see Appendix A) as 77 = y^me/a 2 . The numerical 
values may depend on the specific geometry (triangular lattice vs. square lattice, pore length L, etc.). 

In Fig. [7] we examine the dependence of the translocation time exponent a on the pore diameter. Since the exponent 
a depends on the chain length, we average a(No) over different chain length ranges, giving the average (a(No)) n as 
a function of pore diameter d. For extremely long chains, the size of the pore affects the translocation dynamics only 
slightly. This is because for long chains, the overall friction of the system is dominated by the friction between cis side 
subchain and the solvent, with the pore friction adding only a minor contribution. For shorter chains, however, the 
exponent a clearly decreases with decreasing pore size. For short chains and narrow pores, the pore friction dominates 
the solvent friction, and the dynamics approaches the constant-friction limit, where a = 1 |10j . Similar behavior was 
seen using MD simulations in Ref. [IS], where a was measured as 1.35, 1.30 and 1.21 for 64 < Ao < 256 and pore 
diameters 1.5, 1.3 and 1.1, respectively. Using the same parameter values for the BDTP model, we obtain the values 
1.37, 1.32 and 1.24 for a. The values agree within the statistical error, although there seems to be a systematic error 
of about 0.02. This may be due to a slightly different pore geometry (triangular vs. our circular) used in Ref. [35] . 

To further illustrate the finite size effect of the pore friction, the effective exponent a(Ng) is solved for different 
values of the pore friction r) p . The results are shown in Fig. ^1 The exponent a(No) asymptotically approaches 1 + v 
for all fj p . However, the values for finite Nq and the rate of the approach depend on fj p . For small fj p , i.e., wide pores 
or very viscous solvents, a(N ) is larger and the asymptotic limit is reached for shorter chains than for large rj p (tight 
pores or low solvent viscosities). Similar results were obtained in the MD simulations of Ref. [3B]. The data of Fig. [8] 
indicates that for typical parameters and in the chain length regime relevant for simulations and experiments, the 
pore friction fj p can be a major contribution to the total friction and significantly changes the translocation dynamics. 

C. Dependence of r on the driving force / 

In addition to the exponent a, another important measure of translocation dynamics is the dependence of r on 
the driving force /. Usually, this dependence is characterized by the scaling exponent <5, defined as t ~ The 
simplest argument, namely translocation at constant velocity vq ~ /, gives the scaling exponent S = —1, which is also 
supported by some theoretical studies [U [TUl 021 [20l |2T] , numerical simulations [23 , 26H2H HOI EI] and experiments [6] . 




- ,0 V 



a 



P 



O A 
& 

A' P 

A ei ? 
ff V 



o v 



□n p =i 

Ofi p =2.5 
Afi p =5 

o?i =10 



10 




f 

FIG. 9: (Color online) The dependence of r on / measured from the BDTP model with N = 128 and k B T = 1.2, 7 = 0.7. 
The inset shows the exponent 8 (r ~ /*) as a function of / for both BDTP and MD. Here MD (1) corresponds to 77 = 0.7 and 
k — 15, MD (2) to 77 = 10 and k = 15 and MD (3) to 77 = 10 and k = 150. For S, the error is of the order of the size of the 
symbols. 

For the deterministic tension propagation theory, Sakaue predicts that 8 = — 1 for the SS and equilibrium regimes, 
5 pa —1.26 for the TR regime and 8 pa —0.74 for the SF regime. In Refs. [201151], on the other hand, the exponent 
5 = — 1 is predicted even for the TR and SF regimes. In most MD simulations, the measured exponent is slightly 
larger than —1, typically between 5 sa —0.9 and S ~ —0.97. For the BDTP model, we have measured S pa —0.94 for 
No = 128, ksT = 1.2, 77 = 0.7 and 0.25 < / < 15, which is in good agreement with both experiments and simulations. 

However, looking at 6 more closely reveals interesting details. As shown in Fig. [9J for / < 1, S pa —0.9 and turns 
over to 5 pa — 1 for / 3> 1. Similar behavior for the low-force regime was reported in [34 . However, in the high force 
regime the MD simulations of Ref. [33] give 8 pa —0.8 instead of 8 pa —1. On the other hand, the scaling 8 pa — 1 
given by BDTP for large / is in agreement with the theoretical prediction of Refs. [ToTfTTl H"9Tf2"l"] . To investigate 
the matter more closely, we have performed extensive MD simulations. We have used the same parameters (77 = 0.7, 
kgT = 1.2) as in Ref. [33] to allow direct comparison. In addition, we have studied the effect of solvent viscosity by 
running simulations with 77 = 10.0 and the effect of bond strength by using a FENE spring constant k = 150 (hard 
bonds) and k = 15 (soft bonds). 

The results of the MD simulations are shown in the inset of Fig. [9] The effective exponent <5 as a function of the force 
is found by linear least squares fit from three consecutive [log(r), log(/)] data points. For low to intermediate forces, 
we measure 8 pa —0.9, in agreement with Ref. 34 and the BDTP model. For the large forces, we have —0.8 < 8 < —0.7 
for 77 = 0.7 and k = 15, also in agreement with Ref. [33]. However, for the large driving forces the MD results depend 
on the friction 77 and the spring constant k. For the smaller 77, one has significantly larger 8, and similarly for the 
spring constant k. The reason for the former is that for very low friction, the response of the system to the force is 
not linear due to the inertial term in the equations of motion (see Appendix A). When the friction is increased, the 
mass term becomes less significant and the exponent 8 decreases. Similarly, for extremely high forces and small k, 
the bonds can be significantly stretched by the driving force. This increases the exponent 8 for small fc, as shown in 
Fig. [9] However, we have confirmed that even using overdamped dynamics without the inertial term (see Appendix 
A) and sufficiently hard bonds, the exponent does not reach —1, but stays between— 0.9 < 8 < —0.95. Therefore, 
although the low friction and soft bonds typically used in MD simulations contribute significantly to the difference 
between BDTP and MD, they do not explain it fully. The most probable cause for the remaining difference is the 
absence of the trans side subchain from the BDTP model. For sufficiently high /, significant crowding of monomers 
close to the pore on the trans side occurs, which could increase the scaling exponent 8. For very long chains, the 
effect should be small, because the high friction due to the long tail on the cis side leads to slow translocation even 
for large /, and the friction due to crowding becomes less significant. 

D. Time evolution of the translocation coordinate s 



Finally, we look at the time evolution of the translocation coordinate s(t). In Fig. 10 we compare the results given 
by the BDTP model with MD simulations for a fixed chain length Nq = 100 with different driving forces. Overall, the 
BDTP model agrees very well with the MD results. The curves practically overlap, except for the beginning of the 
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FIG. 10: (Color online) The translocation coordinate 5 as a function of time for the BDTP model (filled symbols) and MD 
simulations (open symbols) for iVo = 100. Other parameters are the same as in Fig. [3] Inset: s(t) for 77 = 0.7 and 77 = 10, 
showing the difference between the inertial and overdamped regimes. 



process, where the s measured from MD simulations lags behind the BDTP solution. This difference increases with 
the driving force /. However, as shown in the inset of Fig. [TUJ the difference disappears when the solvent friction is 
increased from 77 = 0.7 to r\ = 10.0. This shows that the difference is caused by the inertial effects in the low-friction 
MD simulations: here the monomers of mass m need to be accelerated for a time t &cc oc 771/77 before they reach the 
friction-limited velocity. 

The time-dependence of s can also be characterized by the scaling exponent /3, which we define as (s(t)) ~ t& |48j . 
For the chain lengths 100 < No < 500 we obtain f3 = 0.75, which agrees with the MD simulations of Ref. [3T] within 
the statistical accuracy, and is also similar to the value reported in Ref. |28j , where slightly different parameter values 
were used. We have also solved (3 from the BDTP model for much longer chains. The results are shown in Fig. [IT] 
for A*o < 10 6 . Once again, the data shows the extremely slow approach to the asymptotic limit. Even for Nq = 10 , 
the value of j3 continues to decrease, although extremely slowly (not shown). 

In addition to /?, we examine the product a/3. Since at the moment of translocation, s(r) = Nq, the exponents a 
and P are related by a/3 = 1. The relation can be exact o nly in the asymptotic limit N 00, since for finite chain 
lengths, the scaling s(t) ~ t@ is not exact, as shown in Fig. 10 However, we expect the relation to hold approximately 



even for finite Nq. The results are shown in Fig. 11 For short chains, the product quickly increases from a/3 « 1 to 
a/3 w 1.07 at No ss 100, where it attains its maximum value. Thereafter the product slowly decreases, approaching 
the theoretical asymptotic limit a/3 — 1. Also in this case, the approach is extremely slow. For Nq = 10 10 we have 
measured a/3 ps 1.015, with the value still gradually decreasing towards 1. In the finite chain length regime, the BDTP 
model is in excellent agreement with the available MD simulation data. For 100 < Nq < 500 we have on the average 
a/3 = 1.065, which matches exactly with the results of Ref. [31 . 

Finally, we look at the fluctuations of the BDTP model. To do this, we have measured the fluctuations of the 
translocation coordinate (5s 2 (t)) = (s 2 (t)) — (s(t)) 2 as a function of time for both the BDTP model and directly from 
MD simulations. In the BDTP model, the trans side subchain is not included, and therefore we expect that the total 
amount of noise is smaller in the BDTP model than in MD. This is indeed the case, as is shown in Fig. [12} Initially, 
(5s 2 (t)) scales similarly for the BDTP and MD, but as the chain translocates and the length of the chain on the 
trans side increases, the fluctuations in the MD simulation start to increase faster. The effect of the missing trans 
side fluctuations can be also seen in the distribution of the translocation times. As shown in Fig. [12] the standard 
deviation of the translocation times is roughly two times larger in the MD simulations, as compared to the BDTP 
solutions. Including additional noise from the translocated part of the chain may bring both the scaling of (5s 2 (t)) 
and the distribution of translocation times closer to the MD results. However, a proper implementation of the trans 
side fluctuations would also require considering the out-of equilibrium effects of monomer crowding in front of the 
trans side pore entrance. This is a subject of future study and is not within the scope of the present work. 



12 



0.8 


P 


\ 

h 




0.75 




k 




o 
Z 




k 

□ 




CO. 

0.7 








0.65 




\ 









12 3 4 

iog 10 (N ) 




2 3 4 5 

'°9io( N o' 



FIG. 11: (Color online) The exponent /3 ((S(t)) ~ i' 3 , left panel) and the product of the exponents a and /3 (right panel) 
as a function of chain length. The dashed lines indicate the theoretical asymptotic values /3(iVn — > oo) = 1/(1 + v) and 
a/3(No — > oo) — 1. Model parameters are the same as in Fig. 10 
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FIG. 12: (Color online) The fluctuations of the translocation coordinate s (left) and the distribution of translocation times 
(right) for the BDTP model and molecular dynamics simulations. In the BDTP model, the trans side subchain is not modeled, 
which reduces fluctuations. 



E. Initial stages of translocation 



Lastly, we wish to discuss the initial stages of the translocation process. Several authors have suggested that in 
the early stages of translocation, an initial tension blob should be formed at the pore entrance before the tension 
propagation process can begin jTTJ HOI HI] • The theory predicts that the time it takes for the initial blob to form, T; n it , 
decreases with increasing force, as the size of the blob gets smaller. In MD simulations, this process should be visible 
as an initial period of time during which the translocation coordinate s changes only slightly (say, from 1 to 2). 

We have investigated the prediction by running extensive MD simulations and looking at the time evolution of s 
and the dependence of the corresponding scaling exponent /3 on the driving force. The results are shown in Fig. 13 as 
open symbols. At least for / > 0.6, the initial time evolution of s seems to be insensitive to the driving force. Hence, 
we find no indications of the blob initiation process. In addition, we find that the scaling exponent fj, as measured 
using values in the range Nq/2 < s < N$, is also quite insensitive to the driving force for / > 3, and decreases with 
decreasing force for / < 1. This result is in agreement with the BDTP model's result that for low forces the exponent 
a decreases with increasing force and the approximate relation a/3 « 1. However, the result is contrasted by the 
recent study by Dubbeldam et al. |21j where it is found that ft decreases with increasing force, and where the blob 
initiation process is observed to increase the time that the system remains close to the initial value of s ~ 1. 

In an attempt to resolve the contradiction, we have repeated the MD simulations with an additional reflecting 
boundary condition (RBC) that prevents the escape of the first monomer to the cis side. A similar boundary 
condition was also used in the MD simulations of Ref. [21 . The results are shown in Fig. 13 as solid symbols. In this 
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FIG. 13: (Color online) The translocation coordinate 5 as a function of time from MD simulations. The chain is either allowed 
to escape the pore to the cis side (open symbols, same data as in Fig. 

boundary condition (RBC, filled symbols). The results change drastically, if the RBC is used. 



or such an escape is prevented by a reflecting 
Inset: the scaling exponent j3 



({s(t)) ~ t") as a function of the driving force / both without and with the RBC. Parameters are the same as in Fig. 
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case, the results agree with Ref . [3T] , showing a decrease of /3 with increasing force, and correspondingly a decrease in 
the time that the system remains close to s « 1. The behavior of the translocation process in these two cases (with 
or without the RBC) is therefore qualitatively different, especially at low driving forces. 

We argue that the difference is due to a fundamental change in the system's free energy, introduced by imposing the 
reflective boundary condition. For the canonical translocation problem, where the chain is allowed to escape to the 
cis side, the free energy has a maximum but no minimum [51 110) . Preventing the escape to the cis side creates a local 
minimum in the free energy where the system can oscillate, changing the problem into a thermally activated escape 
process, similar to the famous Kramers problem [37]. Therefore, the initial period of slow growth of s is not related 
to the blob initiation, but to the thermal motion of the chain in the free-energy well near s ss 1, which precedes the 
eventual escape across the free-energy barrier to the trans side. In fact, according to our MD simulations with the 
RBC, the time that is spent in evolving from s = 1 to, say, s = 2, increases roughly exponentially with decreasing 
force. This is a clear indication of a thermally activated barrier crossing process. In addition, the influence of the RBC 
disappears for large driving forces where the activation barrier becomes negligible, also in support of our argument. 

Finally, we should acknowledge that while our MD simulations show no signs of the blob initiation process, it 
is possible that the process becomes practically observable only for significantly smaller driving forces (/ -C 1). 
Unfortunately, that regime may be out of reach of current MD simulations. 



IV. CONCLUSIONS 



In this work, we have employed theoretical methods and Langevin dynamics simulations to examine the finite 
chain length effects and the resulting non-universalities in driven polymer translocation. We show that the driven 
translocation problem can be accurately described by the non-equilibrium tension propagation mechanism proposed by 
Sakaue [15) . However, finite chain length effects prevail for extremely long chains. Thus, for chain lengths relevant for 
experiments or numerical simulations, these effects need to be taken into account. Especially, the frictional interaction 
between the pore and the polymer has to be included in the theory, if quantitative accuracy is required. This result 
underlines the significance of understanding the interactions between the polymer and the pore also from the point 
of view of fundamental physics, and is especially important for interpretation of experimental data. 

In addition, we discuss the potential limitations of the present tension propagation theory. Although very accurate 
in the typical experimental and simulation regimes, the absence of the trans side subchain from the theory may result 
in less-than-accurate predictions in the extremely strong force regime (/ > ksTN^) and in the extremely low force 
regime (/ < ksTN^ 11 ). In the first case, the non-equilibrium crowding of monomers on the trans side could increase 
the effective friction. In the other extreme, the fluctuations from the trans side chain, which are neglected in the model, 
may become important. In addition, the tension propagation theory of Refs. [15H171IT5H21) is inherently a deterministic 
theory. However, in the low-force regime the tension propagation process becomes increasingly stochastic, instead of 
deterministic. Exploring the implications of diffusive tension propagation is not within the scope of the present work, 
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FIG. 14: (Color online) Left panel: geometry of the pore used in the 3D MD simulations. The pore is formed by placing 16 
monomers of diameter a equidistantly on a circle of radius (d + a)/2, resulting in a pore of diameter d. Right panel: a side-view 
of the pore and the wall, also showing a typical initial configuration of the chain and the direction of the driving force /. 



but is a subject of future study. 
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Appendix A: Molecular dynamics simulations 



The details of the molecular dynamics simulations that we have used for benchmarking are explained in this 
appendix. In the MD simulations, the polymer chain is modeled as Lennard- Jones particles interconnected by nonlinear 
FENE springs. Excluded volume interaction between monomers is given by the short-range repulsive Lennard- Jones 
potential: 



4e 




e for r < 2 1 / 6 cr 
for r > 2 1 / 6 cr 



(Al) 



Here, r is the distance between monomers, a is the diameter of the monomer and e is the depth of the potential well. 
Consecutive monomers are also connected by FENE springs with 



C/feneW = - l -kRl\n{l-r 2 /Rl), 



(A2) 



where k is the FENE spring constant and Rq is the maximum allowed separation between consecutive monomers. 
For the chain, we use the parameters e = 1, k = 15 and Rq = 2, unless otherwise indicated. The main part of the 

wall is constructed using a repulsive external potential of the Lennard- Jones form U cxt = 4e (^ ) 12 — ) 6 + e for 

x < 2 1 / 6 cr and otherwise. Here x is the coordinate in the direction perpendicular to the wall, with x < signifying 
the cis side and x > the trans side. The neighborhood of the pore is constructed of immobile Lennard- Jones beads 
of size a. All monomer-pore particle pairs have the same short-range repulsive LJ-interaction as described above. We 
have verified that using the simple external potential C/ ox t for the wall gives the same results (within statistical error) 
as using a wall made of monomers in fixed lattice sites, at least as long as the interaction between the wall and the 
polymer is purely repulsive. On the other hand, the geometry of the pore itself may have significant effect on the 
results, as is discussed in the main text. 

Similarly to most of the molecular dynamics simulations in the literature [21] [23l I261 - I281 130"! l3"T| 134"! 135] , we take 
the surrounding solvent into account through frictional and random forces. Thus, each monomer is described by the 
Langevin equation of motion 



-V(£/lj + C/fene + C^cxt) - Wi + Ci, 



(A3) 
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FIG. 15: (Color online) The velocity perpendicular to the wall of individual monomers v x for chain lengths No = 100 (upper 
panel) and iVo = 500 (lower panel) as a function of normalized perpendicular distance x from the wall. The location of the 
tension front is given by R. The different plots corresponds to different instances in time, with the average (■) taken over 10000 
independent runs and a time window At = 1. The solid black line indicates the empirical fitting function of Eq. (Bl) with 
b « 1.9. In the BDTP model, the parameter b is not found by fitting, but fixed by global mass conservation. However, the 
resulting numerical value is comparable to the value found in MD simulations. 



where m is the monomer mass, rj is the friction coefficient, Vj is the monomer velocity, U ex t is the external potential that 
gives the is the external force / = — VJ7 cx t in the pore and Q is the random force with the correlations (d(t)C,j(t')) = 

k R T-- 



ji j8{t — t'), where fcgT is the thermal energy. Typically, we have used the parameter values m = 1, rj = 0.7, 
1.2. The equations of motion are solved with the BBK algorithm |3T] with time step St = 0.005. Initially, the 
polymer chain is placed at the pore entrance as shown in Fig. 
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Eq. ( A3 1 is then solved numerically while keeping 



the first monomer fixed until an uncorrelated initial configuration is generated. After that, the whole chain is allowed 
to evolve according to Eq. (A3 1 until the chain escapes either to the cis or trans side. The latter is recorded as a 
successful translocation event. Typically, we average our data over at least 10 4 such events. 

In analyzing the dependence of translocation time r on the driving force / we also used overdamped Langevin 
dynamics to study the translocation dynamics in the limit of negligible inertial effects. In this case, we neglect the 
inertial term in the equations motion, giving the overdamped Langevin equations 



-V(J7 L j + C^fene + U c: 



Ci 



(A4) 



for each monomer. The conventions are the same as in Eq. (A3). Eqs. (A4) are solved with the explicit Euler 



method [42] with time step St < 0.0001. In this case, we average our data over 10 successful translocation events. 
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Appendix B: Velocity profile 

We have measured the velocity profile v(x,t) of the cis side subchain by running extensive MD simulations in 
3D, with a relatively high driving force / = 10. In Fig. [15] we present the velocity profile for No — 500 obtained 
as an average over 10000 independent translocation events and a short time window of At = 1. From the figure, 
one immediately notices that the curves corresponding to different times t collapse onto a single master curve if the 
perpendicular distance x from the wall is scaled by the location R of the tension front, defined as the distance where 
the average velocity goes to zero and, if the velocity is scaled by the maximum velocity. This master curve can be 
approximated by the expression 



tanh b [x/R + 1 

tanh(fr) ^ < B1 > 

where b is the parameter that controls the sharpness of the profile. In principle, b would be nontrivial function of 
several parameters, such that b = b(v, /, ksT, 77). However, in the BDTP model, the parameter b is not found by fitting 
to MD data, but fixed by enforcing global mass conservation, i.e., requiring internal consistency within the model. 

In addition, it turns out that for the purposes of the BDTP model, the exact form of the velocity profile is not 
crucial. The model seems to be very robust with respect to the different forms of the profile, as long as the number 
of monomers is globally conserved. In fact, we also implemented the model with a piecewise linear velocity profile 
such that v(x, t) = v for x > —R + AR and v(x, t) = ^o^p Here AR is a parameter that controls the shape of the 
velocity profile and is fixed by requiring conservation of mass, similar to the parameter b in Eq. |B1| Typically, the 



difference in numerical results given by the piecewise linear profile and those given by Eq. (Bl ) was comparable to the 



statistical uncertainty, although Eq. (Bl) seems to give a slightly better match with MD simulations. Furthermore, 



in both cases, the asymptotic limit of the exponent a is = 1 + v. 
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